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A Genuinely Multi-Dimensional Upwind Cell- Vertex 
Scheme for the Euler Equations 
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The University of Michigan, Ann Arbor, Michigan 48109 

and Institute for Computational Mechanics in Propulsion 
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Abstract 

A scheme for solving the two-dimensional Euler equa- 
tions is developed. It is based on a new scheme for 
the two-dimensional linear convection equation, and the 
Euler-equation decomposition developed by Hirsch et 
al. [1], The scheme is genuinely two-dimensional. At 
each iteration, the data are locally decomposed into 
four variables, allowing convection in appropriate di- 
rections. This is done via a cell-vertex scheme with a 
downwind- weighted distribution step. The scheme is 
conservative, and third-order accurate in space. The 
derivation and stability analysis of the scheme for the 
convection equation, and the derivation of the extension 
to the Euler equations are given. Preconditioning tech- 
niques based on local values of the convection speeds 
are discussed. The scheme for the Euler equations is 
applied to two channel-flow problems. It is shown to 
converge rapidly to a solution that agrees well with that 
of a third-order upwind solver. 

Introduction 

Much of the understanding of modern upwind schemes 
for the Euler equations has come from designing algo- 
rithms for the one-dimensional linear convection equa- 
tion 


As a consequence of this, problems in two or three di- 
mensions are typically solved in a direction-split man- 
ner, with the upwinding directions normal to the faces 
of the computational cell. This leads to schemes that 
are strongly coupled to the grid on which they are ap- 
plied. Discontinuities that lie along grid lines are repre- 
sented properly when treated in this manner, but ones 
that are oblique to the grid are interpreted incorrectly 
by the built-in “Riemann solver” [2]. This suggests the 
need for designing an upwind-differencing scheme for 
the Euler equations that is truly multi-dimensional, and 
therefore less strongly coupled to the grid. The design 
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of an algorithm of this type should be motivated by the 
two-dimensional linear convection equation 

du du du 

s +c 's +c >5 =0 ' 

and by an understanding of the wave-like character of 
the two-dimensional Euler equations. 

Characteristic information has been used in the past 
to formulate schemes. Moretti's A-scheme [3] and 
the QAZ1D algorithm of Verhoff and O'Neil [4] are 
two examples of non-conservative, characteristic- based 
schemes that use grid-decoupled stencils. Conservative 
schemes that are decoupled from the grid are more rare, 
however. Davis [5] has formulated an upwind scheme 
in which the Riemann problem is not solved normal to 
cell faces, but normal to shock waves. Levy et al. [6] 
have extended this work, including other possible up- 
winding directions. Hirsch et al. [1] have developed a 
method of decomposing the Euler equations into a set 
of convection equations. They have formulated a first- 
order scheme based on this decomposition. Roe [2] has 
developed a different decomposition method, based on 
locally decomposing the data into waves. He has for- 
mulated a first-order scheme that makes use of his de- 
composition. All of these conservative algorithms are 
extremely nonlinear. Differences are not taken in grid- 
contravariant directions, but in directions determined 
by local values of the flow variables. In general, the di- 
rections are actually based on derivatives of flow quan- 
tities. For this reason, these schemes are inherently less 
robust than schemes that use the grid-contravariant di- 
rections. 

Most of the upwind schemes used to date are cell- 
centered schemes. While cell-vertex schemes have ad- 
vantages in terms of accuracy [7], the ones that have 
been developed for the Euler equations thus far are 
based on central differencing [8,9] or on the Lax- 
Wendroff scheme [10,11]. In the central-differencing 
version of a cell-vertex scheme, the residual for the cell 
is distributed equally to the four nodes of the cell. In 
the Lax-Wendroff version, this distribution is altered 
by the higher-order terms, so that the nodes receive 
unequal portions of the residual. This can be general- 
ized so that the nodes receive some weighted fraction 


1 



of the residual, where the weight is determined from 
the stability analysis of the scheme. For a convection 
problem, the weights should be such that the residual 
is “pushed” downwind. For the Euler equations, there 
is the added difficulty of determining what variables 
should be convected, and in what directions. 

For the design of a genuinely multi-dimensional up- 
wind cell-vertex scheme, then, the following compo- 
nents are necessary: 


and are approximated here by the third-order accurate 
one-parameter formulas 




u i,Hi 


UjJ + u»+ U 

2 

-^|{(i + 0«K + i,, + (i (3) 
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1. a cell- vertex scheme with a downwind-biased dis- 
tribution for a scalar convection equation in which 
the grid components of the convection speed are 
known; 

2. a method of locally decomposing the Euler equa- 
tions into a set of convection equations; 

3. an extension of the scalar scheme to a system, such 
that mass, momentum and energy are conserved. 

These components are described below for the case of 
two-dimensions. 


Scheme for the Convection Equation 


The heart of the new scheme for the Euler equations is 
a cell-vertex scheme for a two-dimensional convection 
equation 


du _ du du 
dt Cx dx Cy dy ‘ 


(1) 


This scheme is analyzed below. On a uniform Cartesian 
grid, the residual for the convection equation is given 
by 


~24 {(1 + 0y) UiJ+ 1 + (i “ 9y) u *j) * (4) 

where 6 \ and 6y are the centered second-difference op- 
erators 

6% Uij = -2uij 

Slmj = u<+i,i - 2 UiJ + Uj-,., . 


Using these formulas for the cell-face averages, the 
Fourier footprint of the flux integration for a cell is 


T (At Res) = — ^ v x sin ^-^13 


L . A, 

+i/ y sin -£ 


k- cos *k\ 

2 2 J 

f., 0* 30, \1 

^13 cos — cos — j - J j 

a • 0x ( ■ 3/?y „ . 0y\ 

e y u x sin — Ism - J- - 3 sin y J 

{ . w* o ■ Ml 

V sm — -3sm T j] 


+0 r t/ y sin OL 


where the /?’s are the Fourier variables and the i/s are 
the Courant numbers 


Res i+ iJ+i 


Ax ( U,+lj+ * “*•>+$) 

- A^( U<+ i-' +1 ““ <+ > J ) ’ (2) 


c x Ai 

Ax 

c y && 

Ay ‘ 


where the semi-integer index denotes a average over a 
cell-face, i.e. 

1 / r * +l 

1 /‘ Vj+1 


To update the nodes, the cell-centered residual, given 
in Equation 2, multiplied by At, will be sent to the 
nodes (f+lj), (i+l,j+l) and (tj 4- 1) with 

weights lj sw , u > $e , u > ne , and w nU! respectively (see Fig- 
ure 1). The Fourier footprint of this distribution step 
is given by 


T{Dist) 



4* W sw ) cos 


0x 4“ Py 
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These cell-face averages, to fourth order, may be writ- 
ten as 


4“ + u> se ) COS 


0x-3y 
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If a simple forward-Euler time-stepping scheme is 
used, the net amplification factor for the entire scheme 
is 


G(v Xy Vy, Px, 0 y) = 1 + T (A tRes)f (Dist) (6) 

The appropriate values for the 0’s and the u's remain 
to be determined from the stability analysis. The w’s 
correspond to convection directions, and should there- 
fore be determined by enforcing stability for the long 
waves (l3 Xh 3 y — ► 0). The 0’s control the high-order dif- 
ference terms of the scheme, and should therefore be 
determined by enforcing stability for the short waves 
(At, (3y -+ *)• 

Taking the limit of Equation 6 as i3 x ,f3 y —►0 yields 
the constraint 


fc 2 (Vzftx 4 Vyfly) — 

/ X fix 4 fly , 0 X ~ 0 y 

{u, w - Unt) ^ 4 {^nw - w je ) “ 


where k is some real constant. An added constraint 
(conservation) is that the entire residual must be dis- 
tributed, 


4* u sw 4 u> nw 4 u; se — 1 . 


Also, by symmetry, if the v x and v y are such that con- 
vection is directly towards one node, all of the residual 
is distributed to that node, i.e. 


Wne = 1 if v x = Vy > 0 ; 

u> sw = 1 if i/ x = Vy < 0 ; 

uj nw = 1 if —v x = v y > 0 ; 

u> 5C = 1 if ~v x = v y < 0 . 


Combining these conditions gives the distribution coef- 
ficients 


Une — max 


u> $w = max 


Wn 


= max 


u> se = max 0, 


(°’ (8) 

( 0 - *) .. ) 

V ’ K + t'yl + K - Vy\) 

(Vr-Vy) \ 

Kc-lVl/ ’ 


Wx 4 Vy I + 


( 9 ) 

(10) 


These formulas state that the residual is sent only to 
the nodes that define the downwind face, and is dis- 
tributed in a weighted manner between the two nodes 
on that face. For a plane wave moving in one of the 
coordinate directions, the two downwind weights are 
equal, and the scheme reduces to the standard one- 
dimensional upwind scheme. 


ij + 1 i + 1, j + 1 



A short-wave analysis shows that a necessary condi- 
tion for stability is 

(Wx + Vy | - Wx ~ I'ylMMy + Vy0z) > 0. 

For this to hold for ail values of v x and Vy , the constraint 

0£ _ Vy_ 

Ox ~ Vx 

must be met, so that the 0’s must be given by 

0 r = a f Vy 

Qy — oc v x 

For steady solutions that axe independent of At, the 0’s 
must be independent of At. This leads to the choice 

0 . = -7— — 7 (11) 

max (Wx\i l^yl) 

0y = ;V'i) ( i2) 

max(|t/ r |, Wy\) 

where a is a positive parameter of order one. It is 
interesting to note that this says that the 0’s must be 
downwinded , i.e., in the cell-face average calculation of 
Equations 3 and 4, the 0’s must be chosen so as to 
give more influence to the second-difference about the 
downwind node of the face. 

The only parameter in the scheme that remains to 
be determined is the value of a. Figures 3-5 show the 
effect of different values of a on the amplification fac- 
tor G of the scheme. The maximum amplitude of the 
amplification factor over the high-frequency region (see 
Figure 2) 


G max ( Of ) — - max 


G(f3 X y Py)V Xl Vy,0() 
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High Frequency Domain 
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Figure 2: High-Frequency Region Used to Determine a 
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Figure 3: Effect of a on Stability — 0° Wave 







New Cell- Vertex Scheme 
Stability Boundary 



Figure 6: Stability Boundary for Scheme 


is plotted as a function of a for waves traveling at 0°, 
20° and 45°. Based on these results, the value a = 1.0 
was chosen for the scheme. The stability boundary for 
the scheme, with a = 1.0, is shown in Figure 6. 

The locus of the scheme (i.e. the Fourier footprint 
of F(AtRes)F(Dist)) is shown for the 0°, 20° and 45° 
waves in Figures 7-9. The plots are generated by vary- 
ing f3 x continuously and /3 y discretely, which leads to a 
mesh of points within the continuous footprint of the 
locus. The actual locus of the scheme consists of the 
lines in the plot, and all the space between the lines. 
The circular stablity boundary of forward-Euler time- 
stepping is circumscribed about the loci for reference. 
The loci are very different from those of first-order up- 
wind or central-difference schemes. It is the wave that 
is convected at 45° that is damped the best, while waves 
at 0° (or 90°) are not damped well. This can be seen 
clearly in the contours of the amplification factor for 
the scheme, shown in Figures 10-12. 

Some numerical results for a convected Gaussian on 
a 32 x 32 grid are given in Figures 13-15. In each case 
the Gaussian propagates across the grid virtually undis- 
torted. The onset of a zebra instability can be seen in 
the 0° case, as predicted in the stability analysis. The 
amplitude of oscillations in this case is very small (on 
the order of 10“ 4 ). The convergence history for each of 
the cases is shown in Figures 16-18. The Gaussian con- 
vects at almost one cell per iteration, so that the slope 
of the residual curve changes drastically after approx- 
imately forty iterations. The 45° case, which has the 
best high-frequency damping, converges very quickly, 
while the 0° case converges very slowly. Table 1 shows 
the results of a grid-refinement study confirming the 
third-order accuracy of the scheme. 


u 9 — 1.00 and u y — 0.00 
Wave number locus 



Figure 7: Fourier Footprint of Scheme — 0° wave 


v x = 0.90 and u y = 0.30 
Wave number locus 



Figure 8: Fourier Footprint of Scheme — 20° wave 
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Figure 13: Gaussian Convected at 0° 
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Figure 15: Gaussian Convected at 45° 
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Convected Gaussian 
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Convergence History 
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Figure 14: Gaussian Convected at 20° 


Figure 16: Convergence History for 0° Case 
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v s - 0.90 u y = 0.30 £ = 0.25 u x = 0.80 u y =t 0.80 £ = 0.25 

Convergence History Convergence History 




Figure 17: Convergence History for 20° Case 


Figure 18: Convergence History for 45° Case 


Grid 

Peak at x = 0.5 

Estimated Order 

4 2 

0.84602 

- 

8 2 

0.95929 

1.92 

16 2 

0.99314 

2.57 : 

32 2 

0.99907 

2.88 

Oi 

to 

0.99988 

2.95 


Table 1: Order of Accuracy Study — Convection of 
Unit-Ampilitude Gaussian, £ = 0.25, at 45° 


To solve the Euler equations with a scheme analagous 
to the one above, the system must be decomposed into 
a set of two-dimensional convection equations, with or 
without source terms. Once the equations have been 
decomposed, each component can be treated with the 
convection scheme described above. The distribution 
step carries over in a very straightforward manner; the 
flux calculation (particularly the higher-order terms) 
must be treated carefully to ensure that the formulation 
for the system is consistent with the formulation for 
the scalar equation, and that the resulting scheme is 
conservative. 


Scheme for the Euler Equations 

The Euler equations are 


<9U dF dG 

dt + dx dy 

where U is the state vector of conserved variables 


(13) 


P 


pv 

. ? E > 

and F and G are the flux vectors 


t y 

pu 


r 

pv 

pu 2 + p 


puv 


► G = < 


puv 


pv 2 +p 

pu(E + p/p) ^ 


pv(E + p/p) 


Decomposition of the Euler Equations 


Roe [2,12] has formulated a decomposition of the two- 
dimensional Euler equations, based on the eigenvectors 
of the matrix 


k dF , 3G 

A = au c "'’ + 9u ,,n « 


The eigenvectors of A represent a shear wave, a contact 
discontinuity and two acoustic waves. Roe makes use 
of these eigenvectors, decomposing the flow into, for 
example, four acoustic waves, one shear wave and an 
entropy wave. He uses local values of the flow gradients 
to compute the strength and angle of inclination of each 
of the waves in his model. 

Hirsch et al. [1] have formulated a different decompo- 
sition, which converts the Euler equations to the form 


dW* 

dt 


-f Dx 


dvr 

dx 


+ Dy 


gw* 

dy 


= S, 
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where W* is a vector of convected quantities (entropy, a 
component of velocity, and two acoustic-like variables). 





D x and Dy are diagonal matrices of convection speeds, and 
and S is a source term. Similarly to Roe’s decomposi- 
tion, the decomposition of Hirsch is based on the matrix 


0 


A 


3F a dG . B 

cos 8 4 - — sin 8 

dU dU 


d F 

du‘ 


dG 

: * + du Ky 


where k x and k v are the components of a unit propa- 
gation vector. Diagonalization of A gives 


A = PD(k)P _1 


where 


D(«) 


u • k 0 0 0 

0 u • k 0 0 

0 0 u k + c 0 

0 0 0 u k — c 


and P and P” 1 are the left- and right-eigenvector mat- 
trices. The choice of the unit vector k is based on 
local flow gradients. Hirsh chooses k so as to minimize 
the components of the source term S. To minimize 
all components, one needs two different k vectors; 
for the velocity-component convection and for the 
acoustic-like convection. With two k vectors, 


g | -f (4 1} & - (W 3 * + W?) ^ 

- c (4 2) £-^ 2) £) w? 

-c(4 

These equations are general, holding for any choices of 
and Hirsch shows that, in order to minimze 

th source terms, one needs a that is aligned locally 
with the pressure gradient, and a k ^ that is related to 
the strain-rate tensor. That is, is given by 


K 


M = JZ. 

|Vp| 


and k ^ is computed from the velocity derivatives in 
the following way: if 



dudx 
4 d ^ 


then the propagation angle 


tan 9 = 



2 

L dx 


is calculated; otherwise, the possible propagation angles 
are given by 


AW* = l 


Ap - 

Ky^Ali — k^Av 

• Au — 4? 


— K (2) • Au - ^ 


and 


D x 


u 0 0 0 

0 u 0 0 

0 0 11 + «r 2) C 0 

0 0 0 U — Ky 2 ^C 


tan 0 = 


(fs+s?)±\/(£ + ^) ! ~ 4 * f ^ 
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• (14) 


The value of is then 

k ^ — cos 0i + sin Oj . 

The proper branch for Equation 14 is the one that max- 
imizes the inner product This inner product 

appears in the denominator of entries of the transfor- 
mation matrix P* (described below); the two vectors 
therefore must not be perpendicular. This is ensured 
by taking := if the nominal value of the inner 
product is less than 1/10. 

Hirsch ’s decomposition was chosen for this study be- 
cause the matrix P* is square (4 x 4), as opposed to 
Roe’s decomposition, which yields a 6 x 4 matrix. 


Dy 


v 0 0 0 

0 v 0 0 

0 0 V -j“ K)y~^C 0 

0 0 0 v — «y 2 ^c 


Extension of Convection Scheme 

Just as in the scheme for the convection equation, the 
scheme for the Euler equations is made up of two pri- 
mary steps: 

1. a residual calculation based on a flux integral; 
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2. a residual distribution. 

Each of these steps is somewhat more complicated for 
the system, however. 

For the Euler equations, the residual for a cell is given 
by 


with 




av y 


,(*) 


max (l* / f ) |>l* / n fc) |) 


K-es.+ij+i = -j£jTdy~Gdx] , 


(15) 


where A is the cell area and the integral is taken along 
the cell’s boundary dA. The cells are now quadrilater- 
als, with faces lying along curvilinear coordinate lines 
£ = and T) = tfj. The boundary integral of Equa- 
tion 15 is composed of contributions of fluxes normal 
to cell faces. To extend the approximation of Equa- 
tions 3 and 4 so that they apply to the above residual, 
it is necessary to convert these to flux approximations. 
Equation 4, for instance, when multiplied by c x Ay, be- 
comes an expression for the total flux across the cell- 
face centered at (i, j -f 1/2): 

A C x Uij +C x Uij + 1 A 

c x u iJ + xAy = r —Ay 


The contravariant Courant numbers and and 
the transformation matrix P* are defined further below. 
Note that P* and 9^ are defined per cell, and must be 
averaged over neighboring cells to yied a cell-face value. 
The analog of Equation 3 is 


G; 


*+ 2 J 


a« = 




JllEu 

12 l 


2 

+ 


Gi j + G l+ i j 

2 A * X 


with 


_ f Crli * J a 

12 \ 2 y 


+ 2 ( 0 e ).+a > K F »+ij - F i,i) 

- (G.+ij - G «„j) A ( x] } 

A^X “ ^i+lj **“ 

= yi+ij-yij 
a* = ft+i-fi 


and 


With regard to the Euler equations, this translates di- 
rectly into 


Gd£ = F dy - Gdx . 
The matrix 0$ is given by 


F .\>+i Ar ? = 

+*«♦' A„ y G ‘-i A,, 


( F >J + F * j+1 a G v + G «.;+i A 

12 

{ 2 " y 2 

+ | (0 " 

)«J+£ + 1 — Ffj) ArjU 


— (®»,J + 1 ““ ®i,j) A|jX] } , 


&( = P*diag{^ ) }p*” , (17) 


with 




ai/, 


(*) 


max (l* / €* ) |>l |/ n t) |) 


in which the following notation is used: 


A n x 

&r,y 

— ^ij+1 x i,j 

= Vi j+i - Vij 

on 


An 

— lj + 1 ~ Ij 

with 


denotes the flux normal to a cell-face, 


iJ+i 

F drf 

— F dy — Gdx . 

( A ->y) <+ 

iJ+i 


The contravariant Courant numbers i/^ and are 
related to the wave speeds normal to the cell faces. 
Thus we have, for instance, 

A< , 


scaled such that 


The matrix 0 n replaces the scalar 9 y \ it acts as a scalar 
in each of the convection equations generated by the 
transformation matrix P* 


0„ = P*diag {*(*>} P*“\ 


(16) 


— 2 ^’j+i ” “b ij) 


The Cartesian wave speeds c x and c y are the diagonal 
elements of the matrices D x and Dy introduced previ- 
ously, and are evaluated at (i + 1/2, y + 1/2) by using 
cell-averaged state quantities and the cell-centered k 
values. Note that the factor At in Equation 18 drops 
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out in the definition of 9^ k \ Analagous to Equation 18 
is 


with 


~ 2 ( X *+W “ x *d + x *+i,j+ 1 x »\+ij) 

(A (y) i+ ± J+ i = ^ (ifc+u - y.,> + j/«+ij+i - y<,i+i) • 

The transformation matrix P* is given, in columns, 
by 

r i i 


P* 


(i) 


U 

V 

all vl 
2 

0 


P* 


O) 



P* 


( 3 ) 


P* 


(*) 



JL 

2c 


p 


• +C/c£ 1) ) 

2c«ri) .#c( a ) 

p 

(««<*> 


2c/C^ 


•«( 3 ) +CU-K< & >) 


JL 

2c 



ti/cri) 

• — Cici 1 ^ 


vkM 

• — C/Cy 1 ^ 


kM ■ 

K ^ — CU • /cri)) 


The distribution step requires, in each cell, projec- 
tion of the residual onto the columns of the matrix P*, 
giving weights rW, and multiplication of each of the 
resulting vectors by an appropriate time-step: 


"-♦Wi “ S( r< * )P '"') i+ x Jtj < 20 » 

= (A«‘V*'P-“') itij + r ( 21 ) 


where k varies from one to four. In the above, P*W 
denotes the Ar th column of P* , corresponding to the i th 
wave, and 6 t V^ is the portion of time change in the 
state vector caused by the jfc th wave. Each may 

be divided between the two vertices of the downwind 
cell-face, according to the weights of Equations 7-10, 


with v x and i/ y replaced by i/^ and e.g 


„<*> 


„(*) - 


max 


fa 




(k) 

■<*> _i_ „<*) 




■( fc) - *<*> 


As indicated in Equation 21, it is not necessary to 
take the same value of At for each cell, or even for each 
wave. Spatial conservation, and, therefore, the ability 
to find steady weak solutions, is guaranteed by formu- 
lating the discrete residual on the basis of Equation 15. 
Using different time steps in different cells is a well- 
known technique called “local time-stepping;” using dif- 
ferent time steps for different waves is new, and will be 
called “characteristic time-stepping.” Mathematically 
speaking, the use of a non-constant At is equivalent 
to preconditioning the equations. Local time-stepping 
takes away the stiffness due to spatial variations, and 
may be called “spatial preconditioning;” characteristic 
time-stepping removes the stiffness due to the differ- 
ences among the local wave speeds, and may be called 
“wave-preconditioning.” 

In local time-stepping, one chooses, with some safety 
margin, the largest single time-step value that satisfies 
the stability criterion (see Figure 6) for the local values 
of all pairs i/^ k \i In characteristic time-stepping, 
the time-step for each pair is maximized separately. 
The validity of this practice hinges on the assumption 
that, in the steady state, each residual (Equation 15) 
vanishes separately; for a cell-vertex scheme this, how- 
ever, is not generally true. All one can assume is that 
the sum of all residual components sent to a particular 
vertex vanishes in the steady state. To prevent an im- 
balance among contributions 6 t UW for a particular k , 
arriving in a vertex from different cells, both local and 
characteristic time-step values need to be assigned to 
vertices (i,j) rather than cell-centers (t -j- 1/2, j + 1/2). 

The use of characteristic time steps requires special 
provisions near sonic lines, steady shocks and stagna- 
tion points, i.e. in regions where one of the convec- 
tion speeds vanishes. The linear stability criterion then 
allows of arbitrarily large values of At; in practice, 
however, its value must be constrained by a solution- 
dependent upper limit. How to do this robustly in 
the multi-dimensional case is not yet known; for one- 
dimensional flow, some progress has been reported [13]. 

Boundary conditions are imposed at the walls by en- 
forcing a tangency condition at the vertices on the walls 
and zeroing the mass-flux for the faces on the walls. 
Boundary conditions at inflow and outflow are imposed 
by a non-reflecting condition described by Lindquist 
and Giles [14]. 
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Results 


The scheme described above was used to compute 
steady flows in a two-dimensional channel, with cosine- 
shaped walls yielding a 10% constriction at the throat. 
Two different inflow Mach numbers were taken; 

1. = 1.75; 

2. Moo = 3.50. 

Both cases were run on the 64 x 32 grid shown in Fig- 
ure 19. It is worth mentioning here that the full grid 
was used in the calculations: although the final steady 
state is symmetric about the channel axis, the transient 
states are not, due to the asymmetry of the This 
asymmetry of the Hirsch decomposition is suspected to 
have a negative effect on the convergence to a steady 
state; this remains to be investigated. In comparison, 
the decomposition of Roe preserves flow symmetry. 

The results for the first case are shown in Figures 20- 
30. Figure 20 shows the Mach number contours of the 
steady flow. The compression waves caused by the co- 
sine bump, the coalescence into shock waves, and the 
reflection of the shocks are clearly seen. There are 
some oscillations at the shocks, due to the fact that 
the scheme is third-order accurate everywhere. For 
comparison, results from a grid-biased cell-centered up- 
wind scheme [15] are shown in Figures 21 and 22. This 
scheme is linearly third-order accurate but formally 
third-order accurate only for one-dimensional flow; the 
approximate Riemann solver used for upwinding is Os- 
her’s. Figure 21 shows the results produced without 
limiting of higher-order terms, allowing a fair compar- 
ison with Figure 20. The results are very similar, the 
cell-centered scheme being a little less oscillatory, but 
yielding less defined reflected shocks. The results of 
Figure 22 were rendered monotone by the use of Van 
Albadas limiter [16], leading to a clear loss of resolution 
for all waves, especially the reflected shocks. The re- 
maining plots show the details of the shock-intersection 
region. Figure 23 shows the grid in this region; Fig- 
ure 24 shows the Mach number contours. Comparison 
of the two shows that, although the shock is oblique 
to the grid, it is everywhere captured across two cells. 
The vectors (related to the pressure-gradient) and 

the vectors (related to the strain-rate tensor) are 

shown in Figures 25 and 26. These give rise to the con- 
vection directions shown in Figures 27-29. The first 
shows the convection direction for the shear variable 
and the entropy variable; this is simply the stream di- 
rection. The remaining two show the convection direc- 
tions for the acoustic-like variables. The shocks are ev- 
ident in the acoustic directions. Convergence histories 
for this case are shown in Figure 30. The two different 
convergence rates are for: 



Figure 19: Moo — 1.75 Case — Grid 
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Figure 20: M oo = 1.75 Case — Mach Number Contours 


1. constant At (no preconditioning); 


2. a different At in each cell (spatial preconditioning). 
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Figure 21: = 1.75 Case — Grid-Biased Upwind 
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Figure 23: M «, = 1.75 Case — Grid 
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Figure 22: M <*, = 1.75 Case — Grid-Biased Upwind 
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Figure 29: Mqq = 1.75 Case — Convection Directions 
for Second Acoustic-like Variable 
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Wave-preconditioning was also tried, but the lack of a 
precise control for the in regions where low wave 

speeds occur (inside steady shocks) made this calcu- 
lation actually converge slightly more slowly than the 
calculation with spatial preconditioning only. 

The results for the second case are shown in Fig- 
ures 31-41. Figure 31 shows the Mach number contours 
for the steady flow. The shock pattern is similar to that 
of the first case, but the shocks move at a shallower an- 
gle, so that they do not reflect from the walls before 
reaching the outflow boundary. Due to the strength 
of the shocks in this case, some extra damping was 
necessary to capture them without large oscillations; 
this was provided by smoothing the residuals after each 
time-step with a biharmonic operator. Comparison re- 
sults, non-limited and limited, are shown in Figures 32 
and 33. Again, the new cell- vertex scheme gives vir- 
tually the same results as the non-limited cell-centered 
scheme. Figure 33 shows the substantial loss in resolu- 
tion caused by turning on the limiter. The remaining 
plots again show the details of the shock-intersection re- 
gion. Figure 34 shows the grid in this region; Figure 35 
shows the Mach number contours. The and k ^ 
vectors are shown in Figures 36 and 37, the convection 
directions in Figures 38-40. Convergence histories for 
this case are shown in Figure 41. The three different 
convergence rates are for: 

1. constant At (no preconditioning); 

2. a different At in each cell (spatial preconditioning); 

3. a different At for each convection equation (wave- 
preconditioning). 

Due to the relatively high Mach number, the wave 
speeds do not exhibit much of a spread, so that the 
preconditioning gives only about a 20% gain. In this 
case, the convergence history clearly resembles that of 
the scalar equation: after 150 iterations, the residual 
drops rapidly. 

Conclusions 

A genuinely multi-dimensional upwind-differencing 
cell-vertex scheme for a two-dimensional convection 
equation has been formulated, analyzed and tested. 
It has been extended to the Euler equations, giving a 
third-order accurate conservative scheme. The numer- 
ical results thus far are promising. The four essential 
elements of the Euler scheme are: 


Figure 30: M = 1.75 Case — Effect of Local At 
(CFL=0.5) 
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Figure 31: Moo = 3.50 Case — Mach Number Contours Figure 33: Moo = 3.50 Case — Grid-Biased Upwind 
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Figure 32: = 3.50 Case — Grid-Biased Upwind 

Results (Limiter Off) Figure 34: M <» = 3.50 Case — Grid 
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Figure 35: Moo — 3.50 Case — Mach Number Contours Figure 37: — 3.50 Case k ^ Directions 
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Figure 36: A/oo = 3.50 Case — Directions 
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Figure 39: M = 3.50 Case — Convection Directions 
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Figure 41: M <*, = 3.50 Case — Effect of Local At and 
Characteristic At (CFL=0.5) 
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• the calculation of the convection directions; 

• the residual calculation; 

• the distribution step; 

• the time-step calculation. 

Three of the four warrant further study to produce a 
practical scheme. 

The convection directions chosen here were the ones 
derived by Hirsch et al. [1]. As these are based on 
derivatives of the flow variables, the numerical values 
can be very noisy. One method used to minimize the ef- 
fect of the noisiness of the Hirsch k’s was to freeze them 
after the residual had dropped two orders of magnitude. 
This improved convergence considerably. Experiments 
with other directions (e.g. the streamline direction) 
for the k's showed that they could lead to faster con- 
vergence. Such alternative directions did not work in 
every case, however. Other wave models, such as that 
of Roe [2], remain to be investigated. 

The residual calculation derived here gives third- 
order accuracy everywhere, which is not an advantage 
near shocks. It is not yet clear how to modify the resid- 
ual formulas in order to ensure monotonicity. In addi- 
tion, the highest-order terms in this scheme are strongly 
coupled to the choice of the k' s. These terms turned 
out to be destabilizing in regions where the k’s were 
highly oscillatory. In particular, subsonic cases (not 
shown here) converged only very slowly, and the final 
solutions were not smooth unless a more reliable fourth- 
order term (e.g. a biharmonic term) was added. 

Not all of the underdamped behavior of the scheme 
can be traced to its high nonlinearity. As can be seen 
from Figure 10, the basic convection scheme does not 
damp any combination of a high spatial frequency along 
one coordinate with a low frequency along the other 
coordinate, if the convection is precisely in one of the 
coordinate directions. This lack of damping is caused 
by the vanishing of either 0 X or 6 y . Improvement of 
the convection scheme in this respect requres the intro- 
duction of additional finite-difference terms; these may 
actually be formulated as a smoothing term following 
the distribution step. 

Finally, the time-step calculation, aimed at achieving 
optimal convergence, is far from robust. The technique 
of preconditioning by calculating a value of At for each 
convection equation, at each cell-vertex, can lead to 
large improvements in convergence [13,17]. When one 
of the convection speeds is very small, the potential 
benefit of wave-preconditioning is greatest, but so is the 
danger of taking the time step too large. A satisfactory 
analysis of this remains to be carried out. 
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